Electron Flow in Circular n-p Junctions of Bilayer Graphene 
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We present a theoretical study of electron wave functions in ballistic circular n-p junctions of 
bilayer graphene. Similarly to the case of a circular n-p junction of monolayer graphene, we find 
that (i) the wave functions form caustics inside the circular region, and (ii) the shape of these 
caustics are well described by a geometrical optics model using the concept of a negative refractive 
index. In contrast to the monolayer case, we show that the strong focusing effect is absent in the 
bilayer. We explain these findings in terms of the angular dependence of Klein tunneling at a planar 
n-p junction. 
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I. INTRODUCTION 

The interface of an n-p junction (NPJ) of graphene^"^ 
is fully transparent for electrons approaching it with a 
perpendicular incidence. Electrons approaching the in- 
terface at a finite angle are still transmitted with a high 
probability provided that the transition between the n 
and the p regions is sharp enough.^ As proposed recently 
by Cheianov et al.,^ this high transparency of the inter- 
face offers a way to use the graphene NPJ as an electronic 
lens. The refraction of electron rays in this system fol- 
lows Snell's law with a negative refractive index, which 
is a consequence of the fact that the wave vector and the 
velocity of the valence band quasiparticles in the p re- 
gion are antiparallel. In the case of a point-like source of 
electrons on the n side of the interface, the NPJ provides 
perfect focusing of the emitted electrons on the p side 
if kn = kp, where fc„ (fcp) is the wave number in the n 
(p) region. If kn ^ kp, the sharp focus transforms into 
a smeared focus and a pair of caustics. (For a review on 
the theory and classification of caustics see Ref.lO.) 

According to our earlier theoretical analysis, focusing 
and caustic formation also arises in circular ri-p junctions 
of graphene, where the n (p) region is defined as the area 
outside (inside) a circle. Such a device is found to be 
able to focus an incident parallel beam of electrons into 
a certain spot inside the p region, however the focusing 
is imperfect and caustic formation arises even if kn = kp. 

The interband or Klein tunneling^'^ of carriers in 
bilayer graphene^^'^^ is remarkably different from the 
same process in monolayer graphene. Namely, the bi- 
layer NPJ reflects normally incident electrons with unit 
probability.^ This difference leads to the anticipation 
that the patterns of electron flow in planar or circular 
bilayer graphene n-p junctions are distinct from the pat- 
terns in their monolayer counterparts. 

In this work, we investigate the possibility of control- 
ling the electron flow in bilayer graphene by using a gate- 
defined circular NPJ. We provide an exact solution of 
the effective Schrodinger equation in the presence of a 
step-like circular potential barrier. Using the exact wave 



functions we demonstrate that in contrast to the mono- 
layer case, the focusing of a parallel electron beam is not 
possible in the circular bilayer NPJ. However, we find 
that caustic formation remains a sizeable and possibly 
observable effect even in the bilayer. We also calculate 
the angular dependence of transmission probability in a 
planar NPJ, and use the results of this calculation to in- 
terpret the absence of focusing and the presence of caustic 
formation in circular junctions. 

The paper is organized as follows. In Section II we 
solve the effective Schrodinger equation modelling the 
circular NPJ in bilayer graphene using the method of 
partial waves. In Section III we calculate the angular 
dependence of transmission probability in a planar NPJ, 
and discuss the results of Section II in terms of the trans- 
mission probability function. In Section IV we discuss 
the validity of the model we use, give a brief overview of 
related experiments, and provide a short conclusion. 



II. ELECTRON FLOW IN A CIRCULAR n- 
JUNCTION 




FIG. 1: (Color onhne) An incident plane wave of electrons in 
bilayer graphene is scattered by a circular n-p junction created 
by a gate-induced circular potential barrier V{r). In the n (p) 
region the Fermi energy lies in the conduction (valence) band. 



2 



In this Section wc consider the scattering of an electron 
plane wave on a circular n-p junction in bilayer graphene 
(see Fig.l). Our goal is to calculate the exact scattering 
wave function as a function of system parameters and to 
identify the characteristics of the flow of electrons inside 
the circular p region. We will focus on the case when the 
radius of the circle is much larger then the electron wave- 
length, since in that regime we expect a correspondence 
between results of the quantum mechanical model and a 
simplified description based on principles of geometrical 
optics. 

To model the two-dimensional electron flow in bilayer 
graphene, we use a two-component envelope function 
Hamiltonian which has been derived by McCann and 
Fal'ko.^^ The derivation of this Hamiltonian starts from 
a simple tight-binding model of bilayer graphene, which 
contains only nearest-neighbor intralayer and interlayer 
hopping matrix elements. These hopping matrix ele- 
ments are usually denoted^'^^ by 70 and 71, respectively. 
This tight binding model predicts that the valence and 
conduction bands of bilayer graphene are touching at the 
K and K' points, i. e. at the two nonequivalent cor- 
ners of the hexagonal Brillouin zone. The Fermi energy 
of undoped bilayer graphene lies exactly at the energy 
corresponding to these touching points. In the vicin- 
ity of the K point, the dispersion relations describing 
the conduction and valence bands are both quadratic, 
E±{k) « ±n'^{k-Kf/2m, where m = fg^ w 0.054mo 
is the effective mass. Here the + (— ) sign refers to the 
conduction (valence) band, a is the lattice constant and 
mo is the free electron mass. In this low-energy regime 
the quasiparticle wave functions can be characterized by 
the eigenfunctions of the 2x2 effective Hamiltonian 
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(1) 



where p± = {px ± ipy)- Note that this Hamiltonian de- 
scribes the valence band and conduction band states si- 
multaneously. Similar statements are true for the vicinity 
of the K' point. We note that the dispersion relation and 
the effective Hamiltonian becomes more complex if one 
includes second-nearest-neighbor interlayer hopping ma- 
trix elements in the tight-binding model. In particular, 
such terms lead to a trigonal warping^^"^^ of the quasi- 
particle dispersion. We comment on the significance of 
trigonal warping in Section IV. 

We model the gate-defined circular potential barrier by 
a step-like potential V{r) = V()Q(R — r), where O denotes 
the Heaviside function. Hence the complete Hamiltonian 
of the system under study is 



H = Ho + V{r)l, 



(2) 



where 1 is the 2x2 unit matrix. The validity of this 
model will be discussed in Section IV. Note that the 

same model was used recently to calculate the lifetime of 
quasibound states in a similar system. 



We concentrate on the regime where the potential bar- 
rier forms an n-p jimction, i.e. the Fermi energy Ep of 
the electrons lies between the Dirac point of the bulk and 
the top of the potential barrier (0 < Ep < Vq). In this 
case, the region outside the circle of radius R contains 
electrons in the conduction band (n-typc^), whereas the 
region inside the circle contains holes in the valence band 
(p-type). 

Our aim is to consider the scattering of an incident 
electron plane wave coming from the n region along the 
a;-axis as shown in Fig. 1 and having energy Ep. Such an 
electron has the following two-component wave function^: 



<^(x,2/)=e^'="- — 



1 



1 

-1 



(3) 



where fc„ = \/2mEp /h. In order to derive the wave func- 
tion describing the scattering of this plane wave, we first 
treat the scattering of cylindrical waves and then uti- 
lize the fact that the plane wave (j)(x,y) is a certain lin- 
ear combination of cylindrical waves. Here we note that 
scattering theory has been used recently to predict trans- 
port properties of disordered^''"^'^ and ballistic^^ bilayer 
graphene structures. 

The system has a circular symmetry around the origin, 
hence the Hamiltonian commutes with a pseudo angular 



momentum' operator 



-ihdm + ftcTx, where 9^, is the 



derivative with respect to the angular polar coordinate 
and Uz is the third Pauli matrix. The presence of this 
symmetry simplifies the forthcoming calculations. 

Using the properties of Bessel functions^^ it can be 
shown that in the n region, for any integer j the wave 
functions 



h«(r,^) 



kj(r,v>) = 



Hf\{knr)e-''^ 
K,_^{kr,r)e-^^ 



(4a) 

(4b) 
(4c) 



are simultaneous eigenfunctions of H and Jz , with eigen- 
values Ep and hj, respectively. We denote the ra- 



dial polar coordinate with r. Here Hrr? , Hm' and Km 
denote Hankel functions of first and second kind and 
the modified Bessel function which is bounded for large 
arguments, respectively. There exists a solution similar 
to those in Eq. (4), containing the modified Bessel func- 
tion Im- We disregard it because diverges for large 
arguments. Analysis of the quantum mechanical current 
density in state hj^'' (h^-^'') shows that it is an outgoing 
(incoming) cylindrical wave. On the other hand, kj is an 
evanescent cylindrical wave which does not carry current 
in the radial direction. 



r(2) 
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Inside the circular p region, the regular eigenfunctions 
of the Hamiltonian H having energy Ep are 



/j_i(fcpr)e""^ 



(5a) 
(5b) 



Here kp = y^TTT^Vo — i?F)/ft and j is an arbitrary inte- 
ger. Similarly to the wave functions in the n region, 
and ij are eigenfunctions of with an eigenvalue Hj. We 
disregard other eigenfunctions of H which are divergent 
at the origin. 

Now we consider the scattering of a single incoming 

(2) 

cylindrical wave, h^- . Since [H, J^] = 0, the pseudo an- 
gular momentum does not change during the scattering 
process, therefore the complete wave function describing 
the scattering can be written as 



(1) 



(P) 



(6a) 
(6b) 



in the n and p regions, respectively. The coefficients Sj, 
Aj , Bj and Cj have to be determined from the boundary 
conditions at the interface of the NPJ: the wave functions 
and their derivatives have to be continuous at r = i?. 
Due to the two-component nature of the wavefunctions, 
the two boundary conditions result in an inhomogcncous 
linear system with four equations and the four coefficients 
as unknowns. This system can be solved analytically. 

Having the coefRcients Sj, Aj, Bj and Cj in hand, one 
can determine the wave function describing the scattering 
of the plane wave 4> in Eq. (3) . Making use of the fact 
that22 



gzfcx ^ J2 i"'J^{kr)e' 



(7) 



mG2 



it can be shown that the plane wave can be written as a 
linear combination of incoming and outgoing cylindrical 
waves: 



(8) 



This expansion allows us to use the coefficients deter- 
mined from the analysis of partial waves to derive the 
wave function describing the scattering of the plane wave. 
In the n region. 



and in the p region 



(9) 



(10) 



The complete wave function tp is constructed by tailoring 
and It is built up from cylindrical waves hav- 

ing energy Ep, therefore ip is also an energy eigenstate 
with energy Ep. Since the cylindrical waves fulfill the 
boundary conditions at R, ^ also fulfills them. Finally, 
since in the n region ip contains only the plane wave and 
outgoing and evanescent cylindrical waves (no incoming 
wave), we conclude that -0 is the wave function which 
describes the scattering of the incident plane wave. 




FIG. 2: (Color online) The spatial dependence of the intensity 
of the wave function |'!/;(r)p is plotted in the scattering area. 
Here knR = 300, and kpR — 300 corresponding to n = — 1. 
The solid (dashed) line corresponds to the caustic for p = 1 
{p = 2), where p denotes the number of chords inside the 
NPJ." 

Numerical results for the spatial dependence of the 
magnitude of the complete scattering state (|i/'(t')P) are 
shown in Fig. 2 and Fig. 3 for two different set of param- 
eters. In both cases, well-defined patterns of the electron 
flow can be identified, the wave function magnitude is 
sharply peaked close to the solid curve. This effect is 
almost identical to the one predicted for circular NPJs 
of monolayer graphene."'^^ As we will argue in Section 
III, the geometrical optics model developed in Ref.4 and 
Ref.ll for single layer graphene is also applicable for bi- 
layer with certain restrictions. According to the referred 
theories, the refraction of the incident electrons is gov- 
erned by Snell's law with a negative refractive index. Af- 
ter the refraction, the electrons enter the p region of the 
junction, and the envelope of the electron rays form a 
caustic. These caustics can be identified in the quantum 
mechanical charge density, as it is revealed by Figs. 2 
and 3. 

Despite the apparent similarities of the monolayer and 
bilayer case, the analogy is not complete. In a circular 
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FIG. 3: (Color online) The same as in Fig. 2 with knR = 200 
and kpR = 300 corresponding to n — —1.5. 



monolayer NPJ the charge density is maximal close to 
the meeting point of the two caustic lines, which means 
that the interface between the n and p regions provides 
strong focusing of the incident electrons. This feature 
is missing in our results for the bilayer. In Section III 
we will show that the absence of focusing is connected to 
a general characteristic of interband (Klein) tunneling in 
bilayer graphene. 



III. TRANSMISSION IN A PLANAR n-p 
JUNCTION 



In this section we study the refraction of electron plane 
waves at a planar n-p junction of bilayer graphene. We 
derive the counterpart of Snell's law for this system, and 
calculate how the probability of transmission depends on 
the propagation direction of the incident electron. The 
obtained results will be used to explain our findings for 
the circular NPJ (Section II). 

The studied system consists of a sheet of bilayer 
graphene in the x-y plane which is n-type for a; < 
and p-type for x > 0. The electrostatic potential which 
creates these regions is modelled by a step-like function 
V{x,y) = VbO(x). We consider a conduction electron 
plane wave incident from the n side of the junction. We 
assume that the propagation direction of the plane wave 
is given by the angle a G [— 7r/2, 7r/2], and it has energy 
Ef {0<Ef< Vo). 

To derive the Snell's law for planar n-p junction of bi- 
layer graphene we follow Ref.4. The length of the wave 
vector in the n (p) region is fc„ (fcp), and the length of 
the corresponding group velocity is u„ (vp). (We assume 



that the plane wave is refracted, and do not consider 
the case of total reflection here.) The incident elec- 
tron has the velocity w„ (cos a, sin a) and wave vector 
fc„(cosa, sin a). At the interface this electron is partially 
reflected with velocity w„(— cos a, sin a) and wave vector 
fc„(— cosa, sina). We denote the direction of propaga- 
tion of the refracted wave by (3, hence the velocity of 
the refracted wave is i;p(cos/3, sin/3). Since the refracted 
wave is in the valence band, its velocity is antiparallel 
with its wave vector, and thus the corresponding wave 
vector is A;p(— cos/3, — sin/3). The translational invari- 
ance of the system along the y direction implies that the 
y component of the wave vector must not change during 
the refraction, i.e. fc„sina — — fcpSin/3, which results in 
Snell's law with a negative refractive index n — —kp/kn, 



sin/3 



n < 0. 



(11) 



This form of Snell's law is identical to the one found 
for monolayer graphene NPJs. Consequently, the math- 
ematical formula describing the caustic lines formed by 
the electron rays in a circular bilayer NPJ is also identical 
to the one derived for the monolayer case. This formula 
is given for the monolayer in Eq. (9) of Ref.ll, and it has 
been used to plot the solid curves in Figs. 2 and 3. The 
correspondence between the description of the electron 
flow in terms of quantum mechanics and geometrical op- 
tics is apparent from the figures: the high-density regions 
of the quantum mechanical wave functions are condensed 
in the vicinity of the caustic line. 

We further investigate the refraction of electrons at the 
n-p interface by calculating the probability of transmis- 
sion as the function of the angle of incidence a. The sys- 
tem is modelled by the Hamiltonian H = Hq -\- VoQ{x). 
The wave functions at the n and p regions can be con- 
structed using the results of Ref.6. In the n region, the 
incident, reflected and evanescent modes are given by 



1 



V'ev,n(a;,?/) 



1 

h{a) 



(12a) 
(12b) 
(12c) 



where kny = fcnsina, knx = kn cos a, Kn = 

kn\/l -\- sin^ a and h{a) — (\/ 1 -I- sin^ a — sina)^. In 
the p region, the refracted and the evanescent waves are 



1 



'l/h{P + n) 



(13a) 
. (13b) 
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FIG. 4: (Color online) Angular dependence of transmission 
probability through a planar NPJ of bilayer graphene. The 
corresponding values of the refractive index n are shown in 
the figure. 



a{b) = arcsin-^. From this result one can express (i) the 
propagation direction of the refracted ray using Snell's 
law in Eq. (11), and (ii) the probability of transmission 
(^refraction) as a function of the impact parameter T{b), 
by combining the results for T(a) and a{b). In Fig. 5 
the darkness of the refracted electron rays reflects the 
transmission probability T(b). The figure indicates that 
the electron rays approaching the junction in the close 
vicinity of the optical axis does not have an appreciable 
probability of transmission at the interface, which leads 
to the complete suppression of the focusing effect (cf. Fig. 
5 of Rcf.ll). 




Here kp^ = fcpCOs/3 and Hp — kpy^l + sin^ /3. Note that 
the refraction angle /3 has to be determined from Snell's 
law [Eq. (11)]. 

The wave function describing the reflection-refraction 
process in the n region is ipn = "^inc + '^V'refl + a''/'ev,n, 
whereas in the p region it is tpp = tiprob + &'0ev,p- One 
has to determine the coefficients r, a, t and b from the 
boundary conditions which match the wavefunctions and 
their derivatives at the interface x = 0. The transmission 
probability as a function of the angle of incidence is 



r(a) 



I,/ x|2 ('0rofr|Wa:|'0rofr 
|^(a)| 



h m 



{'ipinc\Vx\i>inc/ 

p- 
P+ 



where 



(14) 
(15) 



is the a;-componcnt of the current operator. Similarly, we 
found that R{a) ^ |r(a)p = 1 - T{a). 

In Fig. 4 we plot T{a) for three different values of the 
refractive index. A characteristic feature of all the three 
cases is the absence of transmission for perpendicular in- 
cidence, T{0) = 0. This behavior has been predicted and 
explained with the chiral nature of the quasiparticles of 
bilayer graphene.^ With the help of Fig. 5 we argue that 
the absence of transmission at perpendicular incidence 
is responsible for the complete suppression of focusing 
in the circular junction we studied in Section II. Fig. 
5 shows several electron rays approaching the circular p 
region and being refracted at the n-p interface. An in- 
coming electron ray can be characterized by its impact 
factor 6, which is the distance between the incoming ray 
and the optical axis (defined as the line containing the 
horizontal diameter of the circular p region, see Fig. 4 
of Ref.ll). Note that the impact factor of rays entering 
the p region is between —R and R. The angle of inci- 
dence, i.e. the angle between the incoming electron ray 
and the local normal vector of the interface at the point 
of incidence, can be calculated from simple trigonometry: 




FIG. 5: (Color online) Refraction of electron rays at the in- 
terface of the n-p junction. The darkness of the refracted rays 
reflects the transmission probability: a white ray corresponds 
to T = 0, and a black ray would correspond to T = 1. The 
dashed curve shows the shape of the caustic derived from 
Snell's law.^^ Only the rays with no internal reflections are 
displayed. Here the refractive index is n = — 1. 

An other remarkable feature of the transmission func- 
tions plotted in Fig. 4 is that transmission is significant 
(> 0.1) in a wide range between perpendicular (a = 0) 
and nearly parallel (|a| ~ 7r/2) incidence. With respect 
to the circular junction shown in Fig. 5, it means that 
electron rays hitting the NPJ further from the optical 
axis have a significant probability to be refracted, and 
hence, to form caustics along the dashed line of Fig. 5. 
This analysis built on the geometrical optics approach 
and the transmission probability calculations for the pla- 
nar junction provides a qualitative explanation of the 
presence of the well-defined wave function patterns pre- 
sented in Figs. 2 and 3. 



IV. DISCUSSION AND SUMMARY 

During the analysis of circular and planar NPJs, we 
modelled the electrostatic potential as a step-like func- 
tion of position which does not couple different valleys. 
This assumption is justified if a <C (i <C A„, Ap, where a 
is the lattice constant, d is the characteristic length de- 
scribing the width of the transition region between the n 
and p sides of the junction and A„, Ap are the de Broglie 
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wavelengths of the considered quasiparticles in the n and 
p regions. The first relation a <^ d ensures the absence 
of intervalley scattering at the interface, and the second 
relation d <C A„, Ap justifies the usage of a step-like po- 
tential in our model. 

Since the geometrical optics model is expected to cap- 
ture the main features of electron flow patterns only when 
the wavelength of the electrons is much shorter then the 
size of the system, the conditions fc„i?, kpR 3> 1 has to be 
fulfilled to observe a pronounced caustic formation effect. 

In the model Hamiltonian in Eq. (2) we neglected 
the trigonal warping term^^'^'* and used an approximate 
effective Hamiltonian which results in a quadratic dis- 
persion relation. In bilayer graphene trigonal warp- 
ing is strong only at very low energies, when Ep < 
71(73/70)^/4 « l.lSmeV. Here 70, 71 and 73 are hop- 
ping matrix elements of the standard tight-binding model 
of bilayer graphene, and we used estimates for them 
from the review of Castro Neto et al.^ For larger energies, 
the dispersion relation is dominantly quadratic up to an 
energy Ep = 7i/2 « 200meV, where it crosses over to a 
mostly linear dispersion. ^^'^^ Therefore our model should 
give a good description of quasiparticles having energies 
between 1.15meV and 200meV. 

To give a numerical example of parameters which fulfill 
the above criteria, we consider a bilayer graphene circu- 
lar NPJ with Fermi energy Ep — lOmeV, gate potential 
Vq = 20meV, transition region width d = lOnm and p 
region radius R = l/itm. Then kn = kp k. 0.12nm~^, 
An = Ap w 53nm, = kpR w 118 and the refractive 
index n = — 1. For these experimental parameters the 
model we used is expected to give a good description of 
electron dynamics in the circular bilayer NPJ, and the 
relation fc„i? = kpR « 118 1 is expected to ensure the 
strong caustic formation effect in this system. 

Finally, we summarize some experimental results 
which support the feasibility of electron optics devices 
in graphene in general. Control of electron flow in a 



ballistic two-dimensional electron system by means of 
gate-defined potential barriers as refractive elements has 
been realized nearly two decades ago.^^ Direct imaging 
of the electron flow in a two-dimensional electron sys- 
tem has also been carried out applying scanning gate 
techniques. ^^'^^ Scanning tunneling microscopy has been 
used to show that oscillations of the local density of 
states around a static impurity can be refocused to a 
remote location. To realize similar experiments making 
use of the negative refraction index in graphene, a trivial 
prerequisite is the ability to fabricate gate-defined tun- 
able NPJs, which has already been reported by several 
groups.^^"^^ 

In conclusion, we have carried out a theoretical anal- 
ysis of electron dynamics in circular NPJs of bilayer 
graphene. We demonstrated that such a system might 
be used to control the flow of electrons, similarly to 
previously realized and proposed electron optics devices. 
We have pointed out similarities and differences between 
electron dynamics in the circular NPJ of bilayer and 
monolayer graphene. In both devices, electron paths 
form caustics inside the circular p region, and the form 
of these caustics can be described with a geometrical op- 
tics model based on the concept of negative refractive 
index. The major difference is that the strong focusing 
of electrons, which is a characteristic of the monolayer 
device, is completely absent in the bilayer. Our find- 
ings are explained in terms of the angular dependence of 
transmission probability at a planar NPJ. 
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